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Abstract. We study the finite element formulation of general boundary conditions for incom¬ 
pressible flow problems. Distinguishing between the contributions from the inviscid and viscid 
parts of the equations, we use Nitsche’s method to develop a discrete weighted weak formulation 
valid for all values of the viscosity parameter, including the limit case of the Euler equations. 
In order to control the discrete kinetic energy, additional consistent terms are introduced. We 
treat the limit case as a (degenerate) system of hyperbolic equations, using a balanced spectral 
decomposition of the flux Jacobian matrix, in analogy with compressible flows. Then, following 
the theory of Friedrich’s systems, the natural characteristic boundary condition is generalized to 
the considered physical boundary conditions. Several numerical experiments, including standard 
benchmarks for viscous flows as well as inviscid flows are presented. 


1. Introduction 

The subject of this article is the finite element formulation of general boundary conditions for 
incompressible flow problems in a bounded domain ft C (d = 2, 3). The velocity field v and 
pressure p are governed by the Navier-Stokes equations 

dv 

(1) p(— + v • Vv) + Vp — pAv = /, divp = 0, 

ot 

together with initial condition p(0) = vo and constants p > 0 and p > 0. For p — 0 we have the 
Euler equations. 

We consider five types of boundary conditions for wall, inflow, outflow, symmetry and 
characteristic conditions, see Table [lj Depending on whether the flow is inviscid or not, the 
boundary conditions change in nature, e.g., no-penetration versus no-slip in the case of a rigid 
wall. Correspondingly, we subdivide the boundary into dQ = r wa n U Ti n U r out U r sym U r c h ar 
with r sym a hyperplane. In what follows, u — (v,p) and B is a symmetric matrix related to the 
negative part of the Jacobian, see below. In contrast to the first four boundary conditions, the 
physical meaning of the characteristic boundary condition is less obvious, since it corresponds 
to an a priori unknown weighting of the different variables, depending on the definition of B. 
It is however the most natural one for a first-order system in the sense of Friedrich, see for 
example hiesj. Note that the outflow boundary condition is often used in order to limit the 
computational domain by introduction of an artificial boundary r out . 

Our approach for developing a discrete weak formulation is outlined as follows. We distinguish 
between the contributions from the inviscid (Euler) and viscid (Stokes) parts of the equations 
and use Nitsche’s method [26] . which has originally been developed for the Poisson problem; 
it has been extended to the Navier-Stokes equations, see for instance U El 12! . In the last 
cited paper the potential of the method to produce a physically meaningful weighting between 

l 


FINITE ELEMENT FORMULATION OF GENERAL BOUNDARY CONDITIONS FOR INCOMPRESSIBLE FLOWS 



/i = 0 

n> 0 

Inwall 

<5 

II 

o 

v = 0, 

Tin 

V = T D 

V = T D , 

Tout 

II 

o 

dv D 

fi— - pn — —p n, 

on Q 

. ov 

T sym 

o 

II 

£ 

v • n — 0, n—~ x n — 0, 
on 

r char 

o 

II 

cP 

53 

1 

(Pbof -B(u-u d )=0. 

Table 1. Considered boundary conditions 


diffusive and convective terms has been clearly demonstrated by comparison with the strong 
implementation of boundary conditions. This idea, which is particularly interesting for high 
Peclet numbers, has then been extended in [3] to turbulent flows by incorporating a wall law 
into the weak formulation. 

In this paper, we use Nitsche’s method to define a weighted weak formulation valid for all 
values of the viscosity parameter, including the limit case of the Euler equations. Our goal being 
the control of the discrete kinetic energy, additional consistent terms are further introduced in 
the discrete formulation. In order to limit the presentation, we focus here on continuous finite 
element spaces. Furthermore, in this paper we only discuss space discretization. 

The analogous treatment for the convection-diffusion equation has been successfully applied in 
the literature, leading to robustness with respect to the diffusion parameter, see for example pj. 
In contrast to the case of the Navier-Stokes equations, the singular limit (the linear transport 
equation) is theoretically well-understood. Additional difficulties which arise in the present 
situation are the variety of boundary conditions and the coupling between velocities and pressure. 
Moreover, the meaning of robustness is not well-understood, since the incompressible Euler 
equations are known to admit very complex solutions. Their mathematical theory is an active 
topic of research, for example the blow-up in three dimensions [21], or the notion of weak solutions 
[23USES], in contrast to the compressible Euler equations, we cannot use entropies as a roadmap 
for the development of numerical methods. We therefore use the kinetic energy as a guideline, 
making sure that the discrete equations do not generate unphysical growth in energy. 

The summary of the article is as follows. Section [2] is devoted to the inviscid equations with 
the characteristic boundary condition. We write the Euler equations as a degenerate first-order 
system and introduce a balanced spectral decomposition of the flux Jacobian in order to define the 
boundary matrix B in Table [l] The term ’balanced’ refers to the fact that the resulting boundary 
condition has the same dimensioning as the equations 0 in the interior of the domain. 

Then in Section [3] we generalize this boundary condition to the other physical conditions of 
Table [l] by letting the data of the characteristic condition depend on the unknowns. For the wall 
condition, such a technique is often employed in compressible flows, using reflection at a solid 
wall. However, it turns out that additional terms should be introduced in order to control the 
kinetic energy. These terms are consistent, except for the outflow condition in case of re-entrant 
flows, where we add an integral which corresponds to a modification of the outflow condition. 
Modifications aimed to increase stability in this case have previously been proposed mi Ena ei 
from a different point of view. 

In Section [4] we add the viscous terms to recover the Navier-Stokes equations. We first in¬ 
troduce the discrete weak formulation for the Stokes equations, based on a generalization of 
Nitsche’s method. Then we present the weak formulation for the Navier-Stokes equations and 
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we briefly discuss the choice of stabilization terms in light of the balanced scaling of the absolute 
value of the Jacobian. Further, for comparison with the proposed method, we present an alter¬ 
native finite element discretization of the Navier-Stokes equations based on strong enforcement 
of the normal velocity in the discrete space. 

Finally, Section [5] presents various numerical experiments involving standard test cases. We 
use the backward facing step problem and the flow around a cylinder to investigate the behavior 
of the outflow boundary condition. The first example also illustrates the necessity to control 
the kinetic energy. Then the Kovasznay flow is used to investigate robustness with respect 
to the viscosity parameter. As examples for inviscid flows, we consider the standing vortex 
problem, the rotational flow given by Fraenkel, and the impact of a jet. Comparisons with the 
alternative discretization are also carried out for the Kovasznay flow and the jet impact problem. 
Although the presented computations were all based on equal-order Q 1 finite elements with 
SUPG stabilization, our theoretical results carry over to other continuous discrete finite element 
spaces. 

1.1. Notation. Let us first introduce some notation. The outward unit normal to dft is denoted 
by n. In 2D, we also use the notation n 1 - = (— ri 2 ,ni) T . We will frequently write v n for v • n 
and — v — v n n. We denote for x G M the positive part by x + := max{x, 0} and the negative 
part by x~ := x — x + . Notice that (— x) + = — x~ and {—x)~ — — x + . Similarly, for a symmetric 
matrix A , we define A~ := RA~R T and \A\ := R\A\R T if A = RAR T with A diagonal and R 
orthogonal. 

We will frequently use the symbols < (and in order to indicate that a quantity is bounded 
above (and below) up to a positive constant independent of the parameters of interest, such as 
physical parameters and discretization parameters. 

Throughout, we let and W h denote finite element spaces of continuous functions constructed 
on simplicial, quadrilateral, or hexahedral meshes of maximal cell width We denote by dx 
the diameter of the cell K. We denote by Uh(t) — (vh(t),Ph(t)) G x W h the space-discrete 
solution and by ^ = ((frhiXh) G Yh x W h a couple of test functions. 


2. Euler equations with characteristic boundary condition 
We define the matrices 


A\u) 


pVil Ji 

.j? o. 


d 

A n (u) = 5>A*(u) 
2=1 


pv n I n 
n T 0 ’ 



0 

0 ’ 


where / is the dx d identity whereas Ji is the d x 1 matrix of elements equal to dij (1 
Then the inviscid part of 0 can be written as a first-oder system in quasi-linear 


<i, j < d). 

form as 


( 2 ) 


71 r V~^ 47 / \ 

M X + ^ A{U) 8^ = ’ 
2=1 



The natural boundary condition associated to 0 is the characteristic boundary condition, 
which we write here as B(u — r d ) = 0 on dfl = r c h a r with B to be defined. 

In the following, we first address the question of non-dimensionalization of the characteristic 
boundary condition and consequently, the choice of B. The simplest choice B = A~ on r c h a r 
does not yield the desired property, as discussed in the next subsection. Then we give a weak 
formulation which satisfies an energy estimate. 










FINITE ELEMENT FORMULATION OF GENERAL BOUNDARY CONDITIONS FOR INCOMPRESSIBLE FLOWS 


2.1. Balancing of the boundary condition. Let us introduce 


(3) 

and define: 


© = 


I 0 

o e 


9> o 


|A| e :=©- 1 |0A0|0- 1 , 

It is important to note that 


A~ e := 0 -1 (0A0) - 0 -1 . 


(4) 


1 -lA\ e + A-*= 1 -A. 


An appropriate absolute value function for Jacobian matrices has been used in the context of 
compressible flows in [1] in order to obtain a proper scaling of eigenvectors. 

Here, the idea is to impose the characteristic boundary condition associated to <§: 

(5) A n (u)~ B (u - vP) = 0 on = r cha r, 

and to choose 6 in such a way that this boundary condition scales as the Euler equations. 

For this purpose, let us begin by stating the spectral decomposition of the matrix QA n (u)Q] 
the proofs of the two following lemmas are given in the Appendix. 


Lemma 1. The symmetric matrix 


&A n (u)Q = 


pv n I On 
0n T 0 


has real eigenvalues and a basis of orthonormal eigenvectors. For d — 2, the eigenvalues are 
given by 


( 6 ) 


x x pv n _L a/ 46> 2 + p 2 n 2 

A = pv n , A p/m := — ±---• 


and the corresponding right eigenvectors by 


r — 


n 




» 2 + AJm 


^p/m^ 


For d — 3, A p / m and r p / m are the same and X — pv n is a double eigenvalue, of corresponding 


eigenvectors 


and 


with {t\, £ 2 , n} cm orthonormal basis 0 /M 3 . 


Lemma 2. Let = (0,x) and — (<//,x'). Then 

(7) 

and 

( 8 ) A n (u)~ e ip ■ ip' = pv~^ ■ ^ 

One also has that 


I A / \i , , 1 I U 1 A '± . XX'+ 26 2 (p n (p' + (x + pv n (p n )(x'+ pv n (p' 

\A n (u)\e1p ■ Ip = p\V n Pn ■ 4>n + - /,„ 9 9 o - 

V 40 + p 2 V 2 


- (X + Am^n) (x' + Am <Pn 


\j 46> 2 + p 2 n 2 

\A n (u)\ e ip ■ ip ~ p\v n \ (</>„) + (p\v n \ + 6)cp 2 n + + g . 























FINITE ELEMENT FORMULATION OF GENERAL BOUNDARY CONDITIONS FOR INCOMPRESSIBLE FLOWS 


Now, let us discuss the scaling of the Euler equations. Let q = p/p, p = //p and 
x = s x x , t — Sft, v = s v v , q — s q q , (p — Sf(p. 

Then, with C and C the differential operators in ft and ft = s x ft , we wish to have that 
(9) C(v,q) = (<p, 0) C(v,q) = (<p, 0). 

Since the first component of C{v, q) is 

^L+v-Wv + Wq= S ^( — ^+v\/v+ Nvh 

at s x \s t s v at s* J 

it follows that 

_ s x _ 2 _ s v 

S v — | Sq — Sy , Sf — 

St s x 

Then, denoting now by B and B the differential operators associated to the boundary condi¬ 
tions on dft and dft , we would like to choose 6 such that the same property ^ holds true for B 
and B. From the expression of the matrix A n {u)~ & in Lemma [ 2 J it follows that the characteristic 
boundary condition © translates into 


AmF - V D ) n + {p~ p U ) = 0, pv n (V ~ V U )„ = 0. 


D 


Dn 


Therefore, the scaling is balanced if A m = A m s v , where A m = pv n — \ 4Q 2 + p 2 v^ /2 according 


to (|6j). This yields 9 — s v 6. Thus, 6 is proportional to p\v n \ and to yjp\p\ at the continuous level. 

Now, taking into account time and space discretization with parameters d t and dx, a natural 
extension is to require similar scaling for discrete solutions corresponding to dt, dx and dt = stdt, 
dx — s x dx respectively. An additional possibility thus appears: 9 proportional to pdx/d t . In 
order to cover all choices, we take 9 as a homogeneous function 4? of degree one: 


0 = \fj (p\v h , n \, y/p\p h I, ■ 


Although 9 is defined locally on each cell K , we do not use a subscript for readability. Note 
that 9 is not necessarily constant on K. 

The obvious choice 9=1 (such that 0 is the identity matrix) corresponds to the standard 
characteristic condition with B = A~ but it does not yield a correct scaling. 

In what follows, we shall mostly use a general 9. However, we will sometimes discuss the 
choice 


( 10 ) 


o 2 = {pvh,nf + c; 


dt 


pd-K 


Remark 1 . For 9 given in (10), we obtain 

( 11 ) 


\A n (u)\ e tp ■ ip ~ p\v n \(p 2 + 0(p 2 n + y. 


We underline that the constants involved in the equivalence (11) are independent of all physical 
and numerical parameters, especially of dx and d t and are therefore valid for arbitrary CFL 
numbers and on locally refined meshes. 
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2.2. Weak formulation. We now propose a weak space-discrete formulation for <§ and ©>• 
For the moment, we neglect the interior stabilization which will be discussed in Section |4j For 
this purpose, we introduce 

a char (uh)(iph) -=\ E [ ‘ V> h~u h ■ A\u h )^fC\ + \ ( \A n (u h )\eu h ■ 

l chai (u h )^h) := [ [ A n (u h )~ e u D -xp h 

Jn JdQ 

and we consider the space-discrete problem: Uh(t ) GV/jX W h, 

( 13 ) P ■ 4>h +a chai (u h )(xp h ) = l chai (u h )(xp h ) Vxph = ((f>h,Xh) e Y h x W h . 

For a smooth solution u = (v,p) to and ([5]), using integration by parts together with the 
property divp = 0, we have consistency: 

(14) pj^ ^ • </> h + a cha »(VU = l char (u)^ h ) 'iiph € v, x W h . 

One immediately gets the energy balance of the formulation. 

Lemma 3. (Energy estimate for a char y ) One has 

a char (y>/i)(V’ft) = 7 y[ \A n (^h)\e^h -iph> 0 , VV’eVft x W h . 
z Jen 

Clearly, any discrete solution uh to 0 satisfies 

(15) f [ ^vl +a chai (u h )(u h ) = f f ■ v h - f A n (u h )~ B vP ■ u h . 

at Jo . 1 Jn Jon 

Since \A n (ph)\e is positive symmetric, © shows energy dissipation through the boundary 
for / = 0 and vP = 0. 

Next we rewrite a char in order to get close to the standard mixed formulation, by removing 
the derivative on the pressure variables. 


Lemma 4. (Reformulation of a cha,r ) Let uh = (vh,Ph) an d ph — (Ph,Xh)- The form a char can 
be written as 

a chai (u h )(ip h ) = ^ ( ( v h ■ V-u/0 -<f>h-Vh- ( v h ■ \7<j) h ) ^4- Jf Xh divv/l “ Ph div ^) 

4” J (^^h^n^h ‘ (fih 4“ Ph^h,n -^-nij^h) h ‘ Iph'j 

Proof. We remark that 

7 ; I (yPh * Ph + Xhdwv h -p h divp h - v h • Vx/i) = / iXhdivvh-phdiv^h)^ [ (phPh,n ~ XhVh : 

z Jn Jn z Jon 

Thanks to Q and A n {uh)uh — (pvh,n v h + Ph n , v h,n), if follows that 

[ \\An{u h )\eU h • p h = f \(pV h , n Vh • <l>h + Ph<l>h,n + XhVh,n) ~ A n (u h )~ & U h - p h 
Jon z Jon z 

and therefore 

[ X \A n { u h)\o u h * Ph T ~ f (PhPh,n ~ Xh v h,n) = [ 7i v h,n v h ‘ Ph T PhPh,n ~ A n (u}f) & Ufo • ph. 

Jon z z Jon Jon z 

□ 




FINITE ELEMENT FORMULATION OF GENERAL BOUNDARY CONDITIONS FOR INCOMPRESSIBLE FLOWS 


Thanks to the reformulation of a char , we can write that 


(16) 


a char (u /i ,)(Yv J ) - U ar (uhX'iph) = Wv h ) • <f>h~Vh ■ (v h ■ V(j)h) 


+ 


/ (Xhdiwh -p h div<t> h ) - / f-<t>h+ / F(u D ,u h ,ip h ) 

JQ J £1 j OCl 


where the boundary contribution J r (vP- ) ph) is defined by: 


(17) 


Hu 5 Ph) : — 2 I A n {up) | ©R/i * Ph T A n {up) * Ph T 2 (PhPh,n Xh^h,n) 

2 A n (lLji)uji • ph + A n (uh) 0 (r UjP) • ph ~\~ ~^(jPhPh, n Xh^h,n) 

= ’ Ph A PhPh,n T A n (v,h ) 0 (r R/i) * ph’ 

3. Euler equations with general boundary conditions 


The relation (16) is the basis for our definition of the weak form for the Euler equations 
endowed with the five types of boundary conditions described in Table [I] For this purpose, we 
write the other boundary conditions in characteristic form, by replacing vP (given on r c h ar ) by 
some r, depending on the type of the boundary condition, on the available data as well as on the 
unknown itself. Then we reformulate accordingly the boundary contribution F(u,u,p) in order 
to see which terms need to be added to the formulation, in order to obtain control over the kinetic 
energy. Next, in order to allow for re-entrant flows, the outflow condition is further modified. 
Finally, we propose a new formulation with boundary stabilization for the Euler equations. In 
order to focus on the treatment of boundary conditions, we do not discuss in this section the 


interior stabilization. We will come back to this topic in Subsection 4.3 


3.1. Reformulation of boundary conditions. Let us note that all boundary conditions on 
dXt can be written under the same form as on r c h ar , that is 

(18) A n {u)~ & (u-u ) = 0 

with u specific to each type of boundary conditions. On r c h ar , we have u = vP. 

In what follows, we choose u according to the available data on the remaining boundaries. 

3.1.1. Wall and symmetry boundary condition. We use the reflection of r, i.e. 

u — u — (2v n n, 0) T . 

Thanks to relation Q from Lemma [2j we have that 

(19) A n (u)~ & (u - u) • p = 


9 \2 

m -V n (x + 


o 2 + M 


Hence, the wall condition v n — 0 is equivalent to (18). 

3.1.2. Inflow boundary condition. Since the velocity is known, we choose u = (v D ,p) T such that 
u — u— (v D — v , 0) T . Then 


A n (u) 0 (u-u)-ip = pv n (v D -v) 


n 

D \ 


bn + e 2 + xi 


A m {y D 


^)n(X T A m Pn 

= ~PVn(v ~ V u ) ■ (f) + a(v - V L ’) n (f) n - P(v - V D )nX 
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where we have put for abbreviation 
(20) a(v) : = pv~ - 


A m 


(pKI - U 402 + falT 


+ A m 




> 0 , 


( 21 ) 


m- = 


A m 


A, 


0 2 + A^ 


> 0. 


The inflow condition v — v u implies (18); the equivalence holds if v n ^ 0 on T- m . 

3.1.3. Outflow boundary condition. Since the pressure is known, we now choose u — {v,p D ) T , 
such that u — u — (0, p D — p) T . Then 

A n (u)~ 0 (u-u) -ij} = Am , (p D - P)(X + A m A) 


0 2 + XI 

1 


(p - p d )x - Pip ~ P U )<t>' 


D \ 


\J 40 2 + p 2 v 2 


with fl introduced in (21). Again, the outflow condition p = p u is equivalent to (18). 

3.2. Reformulation of boundary contributions. We now compute the term in¬ 

troduced in @> on each boundary.This allows us to see what stabilization terms are needed in 
order to obtain positivity of the form when ijj — u. Without loss of generality, we take here 
u D — 0. In what follows, a and fl are those introduced in (20) and (21) respectively. 

3.2.1. Characteristic boundary condition. Since u = vP = 0, obviously 

T{u,u,u) = ^\A n (u)\ou • u 

so J 7 (r,r,r) is non-negative and no additional term is needed. 


3.2.2. Wall and symmetry boundary conditions. It is useful to note first that 


P 


26 2 + p 2 v 2 n 


— Tri \ OL — --, 

2 ' 2y/40 2 + f?V* 

which in view of Lemma [2] yields that 

l|^4 n (u)|e(v,0) • ((f), 0 ) = ^\v n \v ■ (f> + OLV n <f> r , 


It follows that 

( 22 ) 


V n V ■<t> + p<t> n = iUn(w)|©(v,0) • ((f), 0) 


+ (P<t>n ~ XVn) + pV n U -<f>n+ ( pV n ~ Ct)v n <j> n + X^n- 




So finally, we get thanks to (|19|) and to the relation pv n — a — 02 + x 2 


that 


1 A 3 

F(u,u,u) = -\A n (u)\ e (v,0) ■ (w,0) - g 2 A 2 + pV ~ n 


(L) 


2 + (\ - 2A ™ 

1 e 2 + M 


PV n ■ 


_ , | v 2 / 2A 2 \ 

The terms pv~ (lA) and (1 — g 2+ ™2 J V v n are not necessarily positive. Since they are both 
consistent with the boundary condition v n — 0, in order to control them we shall subtract the 
terms pv~v^ • q^ and ^1 — 0 2^2 ^ X v n from the weak formulation. The same approach is used 
for the other boundary conditions. 
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3.2.3. Inflow boundary condition. We obtain thanks to ( [22] ) that 

= ^\A n (u)\e(v,0) • (0,0) + (pc/) n - + (1 - P)x»n, 

so IF(u, u , u) = \\A n (u)\®(v, 0) • (t, 0) + (1 — /3)pv n . The term (1 — /3)x v n being of indefinite sign 
and consistent, it will be subtracted from the formulation. 

3.2.4. Outflow boundary condition. We now get in view of the relation 0 that 

T(u,u,xp) = ^\A n (u)\e(0,p) ■ (0, x) + ^K\v • </> + (1 - p)p<j> n + pv~v </> 

which gives 

F{u,u,u) = ^\A n (u)\e(0,p) • (0 ,p) + | \v n \v 2 + (1 - P)pv n + pv~v 2 . 

The terms of indefinite sign are (1 — /3)p0 n and pv~v • 0. Note that the latter is consistent 
only under the additional hypothesis v~ — 0. 

3.3. ’Energy’ boundary condition on the outflow. In order to avoid the hypothesis v~ — 0 
on Tout, needed for consistency, and to allow thus to treat re-entrant flows, we modify the 
boundary condition as follows: 


(23) 


Spv n v + pn — p D n 


where 5 is a numerical parameter to be determined later. If v n =0 we retrieve the initial outflow 
condition. 

S pv — 

We now choose u = ((1-— —)v n n,p D ) T , such that 

Am 

/ j_ 3p v n D \T 

u u — ( v n --— V n n,p -p) . 

Am 


Then thanks to Lemma [2] we have 

A n (u)~ e (u - u) • 0 = -pv~y n _L • (f) n ± - 


A, 


( Spv n V n +p- p D )(X m 4>n + x) 


6 2 + Xl 


and the new outflow condition (23) is equivalent to ( |18| ) for any 8. 

We assume p D = 0 for the moment and write the boundary contribution with the help of ([1 
as 


r{u,U,i 0 = ^\v n \Vn ■ 4>n + 


7= 9 2 (Spv n V n +p)x+ £ V n Vn4>n + p0 r , 

a/40 2 + p 2 v 2 2 


-/3(Spv n v n +p)<f> n , 

where again (3 is defined by ©• Putting 

Q(v,p, 8) = p (iu+ + (i - 8)v~) v 2 n + ^(8pv~v n +p)p 

and using 


1 . . 1 _ + 

2 KI = -v n -v„, v n = v^ + v n , 


yields 


F(u,u,u) = ^\v n \(v£} +Q(v,p,8) + (8pv n v n +p) ((1 -/3)v n + (id 2 + p 2 v 2 ) 1,2 p-0 V) • 
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The last term of F(u,u,u) is of indefinite sign but is consistent with the modified boundary 
condition, so it can be controlled as previously, by subtracting it from the weak formulation. 
Next, we choose S in order to control the remaining terms. 


Lemma 5. Let 6 = 1 and 6 > p\v n • Then 

Q( 

Proof. Using vf — v~ — \v n \, we have 


Q(V,p, 1) ~ ^J+P\Vn\vl- 


„ n \ _ 2 P 2 + 2 PP V n v n + Op\v n \v 2 n 
VUj Pi U — 2Q 

Thanks to Young’s inequality and to ( v~ ) 2 = —v~\v n \, we can write for any e > 0 that 

2 P 2 + 2 ppv~v n + 9p\v n \vl > (2 - i)p 2 + p\v n \(6 + epv~)v 2 n > (2 - ^)p 2 + (1 - e)0p\v n \v 2 , 

using 0 > p\v n \ as well as \v n \ + v~ > 0. By choosing ^ < e < 1 we get the lower bound. The 
other inequality is obvious, so the announced equivalence holds. □ 


Remark 2. For 6 given in (10) it immediately follows that the outflow terms are controlled 
by p 2 /9 + p\v n \v 2 , which is not the complete energy norm ^\A u (u)\qu • u, since the term Ov^ 
is missing, see Remark [7j However, our approach allows to control the pressure and normal 
velocities, which is not the case in mmm- 

From now on, we take 5 = 1, for which we get (for an arbitrary p D ): 

F{u,u,'ip)-(pv~v n +p-p D ) (l - (3)(f) n + (40 2 + p 2 v 2 )~^x ~ = |KI v-(f>+^{pv~v n +p-p D 

3.4. Weak formulation. Taking into account the previously developed stabilization terms as 
well as the choice of 0 (10) we now define 

a Eu (u h )(ip h ) ■= 2 \(v h ' Wi) ' <t>h ~ Vh ■ ( v h ■ V(j) h )j + (xh div v h -p h div</> h ) 


+ 

+ 

+ 


2 v h,n v h ‘ T Ph^h^n ^n{ u h) &u h 


L 


r c har 

f 

r*wa llUTsymUTin 


^2 T Xh^h,n T Phf^h^n^ 

- in 

(f\vh,n\vh ■ (frh + )j(pVh, n Vh,n + Ph)Xh) > 

l Eu (uh)(iph) ■= [ f-4>h- [ A n (u h )~ e u D ■ ip h + f (-pv^ n v D ■ 4> h + aVn(j) h ,n - XhVn) 
Jn J r char J r in v ’ ' 

+ f v (j}Pxh-p D <t>h,r. 


Remark 3. We have subtracted the positive (and consistent) term — ^ 2^2 on r wa ii U 

rsym; in order to get the same boundary term as on r- m (provided that v D — 0). 


We consider the space-discrete problem: Uh(t) E x W^, 

(24) p ■ <t>h + a Eu {uh)i^h) = l Eu {u h )(ip h ) V^ = (4,x A )6¥ ft xW A . 


)x+P D 4 >n■ 
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Lemma 6. (Consistency of a Eu ) Let u = (v,p) be a smooth solution to ^ with fi — 0 and with 
the boundary conditions given in Table [7J except that the outflow condition is replaced by (23). 
Then u satisfies 

r o 

(25) p j^. (f)h + a E »(^) = l Eu (u)(iP h ) \/ip h = (0 h , X h) 6ViX W h . 

Proof. Since u satisfies the characteristic boundary condition (18) on dLl with u defined in 
Subsection 3.1, it follows from (14) that 

P L £ ' * h + aChar(u)( ^ } = W* G Y h x W„ 

So we only have to check that 

(26) (a Eu -o d “)(«)(^)-(Z E u_ i char )(u)( ^ ) = 0 G X W h . 

Next, we recall from Lemma [4] that 

(a char - l char )(u)(ip h ) = |((v • Vv) • 4>h - v • (v ■ V(f) h (j + J^Xh&vv -pdiv<j> h ) 

+ [ F{u,u,i) h ) - [ f -<j) h 

Jon Jn 


so using the different expressions of and the definitions of a Eu and Z Eu , (26) is equiv¬ 

alent to 




T27j 


IA ^ , n 2A m x 

)Xh v n 


P v n v n '(<Mn + (1 “ 


e 2 + M 


e 2 + xi 


'Un4 > h,r 


I (1 - P)Xh(v - v D ) n - f (pv n v n +p-p D ) Ul-P)<t>h,n + 

J Tin J r ou t V 


^/ 4# 2 + p 2 v\ 0 


-n)Xh)= 0. 


This equality holds true due to the considered boundary conditions. It translates the consistency 
of the stabilization terms. □ 

Lemma 7. (Energy estimate for a Eu ) For all fh = (c/)h,Xh) W h, one has 

a Eu (V7i)(VU = \ [ \A n {^h)\e^h ■ + 7) [ |AiWvi)|e(</>/i,0) • (</>fc,0) 

/nn\ ^^char ^ r wa iiUF S y m UFi n 

(28) ' 'p 


+ 


f (?l^h,nl • Mt + Q(^>h,Xh, 1)) 

^r out 


L out 

,D M3 


For vanishing data f, v B , p D , any discrete solution uh to (24) clearly satisfies 

Jt J Q 2 V * = _ ° EU (MK) < 0. 

Proof. We use the corresponding expressions of F(fi>h->f , h’>f , h) on each boundary. 


□ 


4. Navier-Stokes equations with general boundary conditions 

4.1. Weak formulation of the Stokes equations. Here, we consider the Stokes equations, 
endowed with wall, inflow, outflow and symmetry boundary conditions. We first define the 
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bilinear and linear forms corresponding to the viscous term with weak boundary conditions in 
the sense of Nitsche: 


Asc {u h M := / 
Jn 


= / Vvfe : V<f> h 


L 


(29) 


L 


^wallCTi] 


dv h 

• 9h + Vh 
on 


d<t>h 

dn 


UK 




\ dn 
dvh 
r out 0 V dn 


H < t > h,n “ 1 “ I'h.n 

d(f>h 


dn 


[ \ ( • n Xh + • np h - p— • n —- • n 

JTnu t U 


dn 


/r 0 ut 


dK n ' n 
dvh dj>h 

dn dn 


1 d(j>h d 
*np , 


6 dn 


where 7 is a stabilization parameter; the terms on wall, inflow, and symmetry boundaries are 
standard, see for instance |3 Q 21 E|, whereas the outflow is treated in order to fit with the previous 
formulation for the inviscid case. 

Then the space-discrete Stokes problem reads: Uh(t) GV/jX W^, 


(30) 

where 
st 


P j Pfr • 4>h + a st (u h ,ip h ) = l st {?ph) vy> h eVfcXWj 


r out 


(Uh, VU := pa vlsc (u h , ip h ) + / (Xh div v h - p h div <f> h ) + / (Ph^h,n ~ XhVh,n ) + / 

Jo, J r wa ^uri n ur s y m Jt* 

l S \^h) ■= pl V1SC (<t>h) + f f ■ (f>h~ f XhVn f P D <t ) h,n+ f Xh■ 

*/ O Tin ^/ Tout U Tout 

The well-posedness of this discrete Stokes problem for 7 sufficiently large follows from standard 
arguments. 

d v 

Remark 4. Here we have not considered the characteristic condition (p—~ , 0) T — A n (u)~ @ (u — 

dn 

u D ) = 0 since it is not natural for the Stokes equations. For the generalization to the Navier- 
Stokes formulation, this condition leads to the following additional terms to a St and Z St respec¬ 
tively: 

/ (Ph4*h,n — An( u h) &u h ' 'Ph) 5 — j ^-n{ u h) 

** I- r* Vi p\ r ** t nVi p\ v 


It is obvious that the formulation (30) is consistent in the sense that a sufficiently smooth 
solution u — (v,p) to the Stokes equations satisfies the discrete equations 

dv 


■/ 

J Q 


p l^-fr'<f>h + a (u, i> h ) = r'i'iph) Viph, = (4>h, Xh) eV k x w h . 


Depending on the employed discrete spaces, one may need to introduce additional stabilization 


terms; we discuss it in subsection 4.3 


4.2. Weak formulation of the Navier-Stokes equations. The scaling of the Euler equations 
has been discussed in Subsection 2 T We are now dealing with the Navier-Stokes equations, so 
we also have to consider the viscous term. Since 


Av 


rAp, 


1 

nPhXh 
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it follows according to © that s v s x = 1. Together with the previous relation s v = s x /st , this 
leads to the well-known unique scaling of the Navier-Stokes equations which satisfies ©, 

2 1 1 1 

S t — s x , S v — S q — ~21 s f — “3 • 

^x v x 

Recalling that 6 = s v 6 , we deduce that 6 may now also depend on /jl/cIk- 
We now take into account both the viscid and the inviscid parts and define: 

cF aSt (u h )(ip h ) :=a Ell k)(4) + j«a VISC (%,4). 

Z NaSt K)(W 0 :=l Eu (u h )^ h ) + ^r sc (cf> h ). 

In order to take into account stabilization, we define 
a Sab t ( u ^)(V’/i) : = a NaSt (u h )(lp h ) + a stab (u h )(lp h ), ■ = Z NaSt (u /l )(V’/ l )+Z st ab(w/ l )(V’/i) 

with the forms a sta b and 4tab to be defined in Subsection |4.3| 

The space-discrete problem then reads: Uh(t) GV/jX W^, 

( 31 ) P ■<t ) h + af t ll t (u h )('ip h ) =lstlb( u h){^h) V4 = (4,x/ 1 )eV) l xWj. 


4.3. Balanced SUPG stabilization. Following the idea of SUPG jTOj 22], we consider the 
following stabilization for the Navier-Stokes equations: 

astab (“ft) WvO := Y ( / lK,\E Uh (u h ) ■ E Uh (ip h ) + / 7^,2 div v h div <j> h ) , 

A-eK: h Jk X 


4tab {Uh)(lph) •= E f lK,lf ' E Uh (ip h ) 

K&tC h JK 


with the parameters 7 K,i to be specified and with 

E Uh (iph) ■= P + pvh ■ ^(Ph ~ pA(p h + 

Clearly, a stah (u h )(u h ) > 0 and a stah (u)(ip h ) = 4tab(«)(U) for u sufficiently smooth solution of 
the Navier-Stokes system. 

In [ 6 ], we have discussed how to tune the stabilisation parameters in order to get robustness 
with respect to the (local) Peclet number Pe = p\vh\dx/P- For large Peclet numbers, the flow 
is governed by the Euler part of the equations, whereas for small Peclet numbers, the Stokes 
part becomes dominant. The parameters 7 x,i depend on the (local) parameter 6 such that the 
stabilisation has the same scaling as the equations. Following [ 6 ], we take next on each cell K 


(32) « 2 ,= (pv k f + + 4 (^) 2 

with the constants cat, cst satisfying c 2 dt + c| t > 0. The stabilization parameters are taken as 
follows: 

Ik , 1 = Ik ,2 = 72 dxO , 715 72 > 0 . 


Note that 9 is not necessarily constant on K. This cell-wise definition induced by the SUPG 
stabilisation is compatible with the one previously introduced in ( 10 ) for the boundary conditions 
in the case /1 = 0 . 

Also note that we thus recover formulas similar to those proposed in the SUPG literature 
(see for instance m). One retrieves the stabilisation used in the inviscid case by simply taking 


(i — 0 . 






F4NITE ELEMENT FORMULATION OF GENERAL BOUNDARY CONDITIONS FOR INCOMPRESSIBLE FLOWS 


Remark 5. Recalling that min{a, 5, c} can be approximated by (1/a 2 + l/b 2 + 1/c 2 ) -1 / 2 , we can 
write that 

. r dx d t d\ rP\ V h\ P P ^ 

7k ,1 ^ 7! mm{——r, —, —}, 7^,2 - 72max{ ——, —, 

P\Vh\ P P d K d t d z K 

4.4. Kinetic energy estimate. By putting together the previous results, we immediately get: 
Theorem 1. (Coercivity o/a NaSt / For all — ((/)h,Xh) £ x W h, one has 


' r wa nur i n 


(iphK'iph) =P [ |V ^| 2 + /x [ 

Jn Jr 

\[ \Ai(^h)\eiph -i>h + \ f 

z J r„w 1 Jr 


7 ,2 0 A d< ^ h 

- 2*. ■ 


+ M 


/, 


7 .2 r,A d ^ h 
d Kn 2 ^h,n Qn 


n 


+ x 


^wa llUr sym Ur in 


|AiW7i)|e(</>/i,0) ' O/i,0) 


+ 


fS0h.nl (<Ph)t ' (0fe)4 + Q(4,Xft ~ P~izr -n , 1) ) • 


Tout 


dn 


Remark 6. For 9 given in (3^ ) and 7 sufficiently large we obtain, according to Remark ^ 

a NaSt 0h)0/i) ZP [ |V0h| 2 + [ FJh+ [ f~4>h,n + \ [ \An(^h)\e^h ■ i>h 

Jn J r wall ur in a K J r sym a K * J r char 

+ 1 f OOmOI + Hh,n) + l [ \p\<t>h,n\Jh + -(xh- ) • 

1 j r wa iiur sym ur in 1 J r out \ V V dn J J 

Theorem 2. (Energy estimate for Navier-Stokes) Any discrete solution uh to (31) satisfies 
~JT [ ~ ~ P [ 4" f f * Vh — ( a stab ( u h)( u h) ~ ^stab{ u h){ u h)) 

at Jn z Jn Jn 

I (\An{^h) lo^h T A n (Uh) ) Ufa 

^ Eliar 

f (, Dn ( 7 dv h \ dv h \ f ( 7 < 9 ^ 

■" L,ur» U -”> ■ ■- U + Vh - U - ■ _ ' 

~\ [ \A n (u h )\ e (v h ,0) ■ (v h ,0) - I pv^ n v D -v h + v^(p h -av hjn ) 

-T wa llUr s y m Uri n '"'Pin 

gK"! iVh)n ■ (A)n + Q( v h,Ph — ■ n > 1)^ + Jr (^ (^ — 


• n — 


For vanishing data f, v u , p u one gets 


(j, I ^1 = —0 NaSt (u/ l )('U/j) - Ostabluft)!^) 
J 


d f 

so — < 0 for 6 given in (32) and for 7 sufficiently large, 

dt Jn 2 


4.5. Conservation of momentum. In addition to kinetic energy, further physical quantities 
such as linear momentum, vorticity, and helicity are conserved by solutions to the Euler equations. 
A natural extension would then be to construct approximations with similar behavior. These 
questions have been recently discussed in [15, JL6J, where for example momentum conservation 
in simple domains has been shown for isogeometric B-spline methods, with the key ingredient 
of discrete solenoidal functions. The purpose of this subsection is to give a discrete balance of 
linear momentum for our method based on classical Q 1 x Q 1 finite elements. 
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The continuous equations have the global momentum balance: 

(33) p^- [ v = [ f T [ q(u ) with q{u) := (pS7v — pi + pv ® v)n. 

dt Jn Jn Jon 

In order to obtain a discrete analogue of (33) we integrate by parts the term pv^ • (vh * V</>/*)/2 
coming from the anti-symmetrization in a Eu . We then have 

(34) a Eu («fe)(V’/i) = ao(uh)(iph) + b(ph,<t>h) ~ KXh,Vh), 
where: 

(35) 


ao(uh)(iph) = [ P\(vh-'Vv h )-(f> h + l-divv h {v h -(f> h )]-[ pv. v h ■ 4> h + [ \{pv h , 

Jn \ 4 / Jaa ’ ,/rout y v 

4“ j aVfi^ n <f>h,n 4“ j \P^h d" Phfyh^n -^-ni^h) ® nh 

Jr wa iiur S y m ur in J f,.i, « T 

b(ph,<f>h) = ~ / Pft div (f>h + / 

Jn Jr 


,Vh,n +Ph Xh 


' r c h ar 
Ph^h.n- 


^ r wa iiur in ur S y m 

Let e z denote the unit vector of the i-th coordinate. Taking the test function i/j l h = (e z , (pp^ • 
e z )/2), which we are allowed to do thanks to the weak formulation and to the equal-order ap¬ 
proximation, we obtain from ([351) for the Euler part and from (29) for the viscous contribution: 


NaSt 


( u h)(^h) -o NaSt (u h )(^) = [ /*+ [ q l (u h )-f e l h {u h ), l<i<d 

./O J dQ J dQ 


with 


- r d ) on T char 

“ P V h,n V f l -(.Ph-p D )n+ (pVh, n v h ,n + Ph ~ P D ~ • nj v h on r out 

£h(«h) := < v h,n (an - %v h ) - *wj" n Uft + %v h 

v h,n (an — %Vh) — pv hn v h + j^Vh, n n — /i ("5^)^ 

, K,n - W°) (a» ~2 V h)~ P v h,n( v h ~ V°) + “ « D ) 


on r wa n 
on T S y m 
on T in 


and with 4/^ the matrix of lines ip l h . Therefore with q h := q — £ h and taking into account SUPG 
stabilisation, we have the discrete momentum balance: 




KeK h ' 


lK lP l P ~dt + PVk ’ V Vh ~ pAVh + V Ph ) ‘ V Vh + f dQ <lh(uh)- 


4.6. Alternative discretization with strong enforcement of normal velocity. In order 
to ease the comparison with other formulations of the Navier-Stokes equations, we integrate by 
parts in a Eu the term *0^/2 and use the expression of A n (uh)u} l -'iph on Tchar to obtain: 

(37) 

ao(uh)(^h) = ~ P (v h ■ (vh ■ V0 fe ) + i div v h {v h ■ (f) h )^j + pv\ n v h ■ cf> h + i [p v h, n v h,n + Phj Xh 

T / P V h _ Xh v h,n T A n ( u h)~^ e u h 

JTnhBT V 


' r wa n ur S ymur i n 


Starting from (37) we give an alternative formulation with strong enforcement of the normal 
velocity on the inflow and wall boundaries, following the ideas of [a m- On the outflow, we 
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impose the ’energy’ boundary condition pv~v Apn = p°n. In addition, in order to highlight the 
influence of the scaling, we also consider a characteristic boundary with the standard condition 
A~(u){u — vP) = 0 with 9 = 1 (that is, not correctly scaled). 

Let the space of the velocity test-functions = { vh E V^; v^ n = 0 on dQ \ (r out U r c h ar )}- 
The alternative discrete formulation is obtained by replacing ao in (37) by 


ao(u h )(lp h ) = - p(l> h -(v h -'V<l>h)+ P v hn V h- ( l ) h+ (~XhVh,n +A+(u h )u h - 1 p h ) , 

J q J r c Vi ar 


Using (34) we obtain the corresponding expression of a^ u and consequently of a 1Na . The viscous 
term and the SUPG-stabilization are unchanged. 

The alternative discrete problem which will be further used for comparison reads: Uh(t) E 

v h X W h , 

(38) p ^ • 4> h + (a NaSt + a sta b)(uh)('iph) = (P aSt + 4tab )(u h )(ip h ) Viph eY h x W h . 


5. Numerical experiments 


First we use an analytical solution of the Navier-Stokes equations for all p > 0 in order to 
validate our code. We next carry out some experiments concerning the outflow condition. In the 
case of re-entrant flow, comparisons with the ’do-nothing’ condition are also made. Finally, we 
present three test cases for inviscid flows. 

We use quadrilateral meshes and Qi x Q i finite elements for the spatial discretization. For the 
time-dependent problems, we use the BDF2 scheme with small time-steps in order to make sure 
that the spatial errors dominate. The values that we use for the different stabilization constants 
are the following: 

c dt — 0.1, c st = 4, 7 = 100, 7 i = 0.25, 72 = 0.1. 

We have discussed these values in 0. 


5.1. Kovasznay’s exact solution. 

5.1.1. Convergence under mesh refinement. In order to validate our code, we use the analytical 
solution given by Kovasznay for a solution of the Navier-Stokes equations, mimicking the two- 
dimensional flow behind a grid of circular cylinders [24]. The solution is given by 


(39) v = (1 — e Xxi cos(27TT2), e Xxi sin(27TX2)) T , P = 


— e 


2\x± 


A = 


87r 2 /x 


1 + y/l + 167T 2 /i 2 


We use the computational domain Q =] — 0.5,10[x] — 0.5,1.5[ and prescribe the analytical 
expression (39) on the left (x\ — —0.5), the upper (x 2 — 1.5) and lower (x 2 — —0.5) boundaries, 
whereas the right (x\ = 10) boundary is treated as an outflow. Three test cases with p = 
0.025,0.0025,0.00025 corresponding to Reynolds numbers 40,400, and 4000 are considered. The 
vector fields is sh own in Figure [l] 

In Tables 2p and [ 4 ] we show the pressure error in L 2 (fl)-norm and the velocity errors in 
and L z {Ct)~ norms; for each one, we also indicate the convergence order 0{h a ). The 
velocity errors in i7 1 (fl) show the expected first-order convergence with slightly increasing values 
for higher Reynolds numbers. Second-order convergence of velocities in L 2 (Q) is observed. The 
velocity errors vary barely with respect to the Reynolds numbers. The convergence in pressure 
is better then first-order, as typically observed for equal-order stabilized methods. It shows a 
surprising decrease for fi = 0.00025, which is probably due to the special form of the pressure. 
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(a) ii — 0.025. 





(b) fj, = 0.0025. 







(c) n = 0.00025. 



(d) /i = 0.0001 at t ~ 30s. 


Figure 1. Kovasznay flow with N — 3072 nodes: velocity fields 
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Further decreasing the viscosity leads to bifurcation. A typical periodic solution for fi = 0.0001 
is shown in Figure [lj 


N 

\\P-Ph II 

order a 

!|V(r-r /( )|| 

order a 

\\v - V h \\ 

order a 

48 

6.38e-01 

- 

5.89e+00 

- 

7.76e-01 

- 

192 

1.34e-01 

2.24 

4.00e+00 

1.31 

2.58e-01 

1.58 

768 

2.83e-02 

2.24 

1.93e+00 

1.04 

5.43e-02 

2.25 

3072 

8.15e-03 

1.79 

9.61e-01 

1.00 

1.39e-02 

1.96 

12288 

2.27e-03 

1.84 

4.80e-01 

1.00 

3.59e-03 

1.94 

49152 

6.09e-04 
Table I 

1.89 

2. Pressu 

2.40e-01 
re and velocity 

1.00 
r errors, / 

9.06e-04 
t = 0.025. 

1.98 


N 

l\p-ph\\ 

order a 

V(n — Vh)\\ 

order a 

\\ v ~ v h || 

order a 

48 

4.13e-01 

- 

7.22e+00 

- 

5.42e-01 

- 

192 

1.64e-01 

1.32 

6.48e+00 

0.15 

4.01e-01 

0.43 

768 

2.03e-02 

3.06 

3.21e+00 

1.01 

8.86e-02 

2.17 

3072 

3.67e-03 

2.46 

1.60e+00 

1.00 

2.21e-02 

2.00 

12288 

1.12e-03 

1.70 

8.00e-01 

1.00 

5.69e-03 

1.95 

49152 

4.14e-04 

1.43 

4.00e-01 

1.00 

1.47e-03 

1.95 


Table 

5. Pressure and velocity errors, fi — 0.0025 



N 

Up- Phil 

order a 

1 V(n — Vh)\\ 

order a 

\\ V ~ V h || 

order a 

48 

1.76e-01 

- 

9.73e+00 

- 

4.76e-01 

- 

192 

6.40e-02 

1.45 

9.22e+00 

0.08 

4.44e-01 

0.09 

768 

1.27e-02 

2.32 

4.50e+00 

1.03 

8.95e-02 

2.31 

3072 

2.66e-03 

2.26 

2.22e+00 

1.02 

2.07e-02 

2.11 

12288 

4.10e-04 

2.69 

l.lle+00 

1.00 

5.20e-03 

1.99 

49152 

1.92e-05 

4.43 

5.52e-01 

1.00 

1.37e-03 

1.92 


Table 4 

. Pressure and velocity 

errors, /j, 

= 0.00025 



5.1.2. Comparison with the alternative discretization. In Figures [2] and [3] we compare the errors 
with the alternative discretization with strongly imposed normal velocity given in (38) for fi — 
0.025 and fi — 0.00025, respectively. The labels ’strong’ and ’weak’ in the figures correspond to 
the implementation type of the boundary condition. We find that the velocity errors are globally 
very similar and that the pressure errors are slightly smaller for the discretization proposed in 
this article. 


5.2. Outflow boundary condition. In this subsection we investigate the behavior of the pro¬ 
posed outflow boundary condition. In practice, outflow boundary conditions are often used in 
order to limit the computational domain. The so-called ’do-nothing’ condition pn — (id n v — 0, 
resulting from the standard weak formulation of the Stokes equations with vector laplacian, is 
widely used in practice, since it is easy to implement and yields satisfactory results in many 
situations. 
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(a) Pressure error in L 2 -norm (b) Velocity error in if 1 -seminorm (c) Velocity error in L 2 -norm 

Figure 2. Comparison of discretization errors for /i = 0.025 



(a) Pressure error in L 2 -norm (b) Velocity error in if 1 -seminorm (c) Velocity error in L 2 -norm 


Figure 3. Comparison of discretization errors for p — 0.00025 


In the case of a cylindrical domain cut by a plane section, it amounts to prescribing the mean 
pressure if no-slip boundary conditions are given on the walls of the cylinder [20]. An explication 
for its success is the fact that the Poiseuille flow in a channel satisfies the boundary conditions, 
and the channel could therefore be cut at any position for this flow. Therefore, the general advise 
is to use a sufficiently long computational domain in order to cut the domain at a point where 
Poiseuille flow can be expected. However, in some applications it may be difficult to determine 
such a distance in advance, or even impossible due to changes in time or changes in parameters, as 
in optimization problems. Therefore, it is important to understand what happens if the domain 
is not cut at a stationary Poiseuille-type profile. 

One major difficulty related to the outflow condition is its stability. The ’do-nothing’ con¬ 
dition has been reported in the literature to lead to instabilities, if the flow is re-entering the 
computational domain [5]. As shown in Theorem [2j we have control over the kinetic energy for 
our formulation. 

Since our proposal leads to a modification of the boundary condition, the purpose of the 
following experiments is to investigate the behavior of the formulation, especially the changes 
induced by the additional terms. 

5.2.1. Backward facing step. In order to illustrate the possible instability of the ’do-nothing’ 
boundary condition, we consider the following non stationary problem. The initial condition 
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is defined by the stationary solution to the backward facing step problem at Re = 800. The 
domain is chosen to cut the second recirculation zone, such that there is re-entrant flow. The 
non stationary computation is done with a ten-times smaller viscosity, such that Re = 8000. We 
use the implicit Euler scheme with small time steps. It turns out that on coarse meshes, the 
’do-nothing’ boundary condition becomes unstable after a certain time. As shown in Figure |4j 
the velocity field blows up at the re-entrant part of the outlet boundary. Consequently, for 
this condition the control over the kinetic energy is lost, see Figure [5| On the contrary, the 
’energy’ boundary condition allows to control the kinetic energy cf. Figure [5j and this even for 
a longer-term simulation. 





Figure 4. Comparison of velocity fields for the ’do-nothing’ condition (upper 
image) and the ’energy’ boundary condition (lower image) at three different times 
close to the blow-up in the first case 

5.2.2. Stationary cylinder benchmark. Our next test is the computation of the drag coefficient in 
a stationary flow around a cylinder in a channel. The geometry, with a slightly non-symmetric 
cylinder, is taken from m However, we have chosen a two-times smaller viscosity (which yields 
a Reynolds number of 40), in order to enlarge the recirculation zone behind the cylinder with 
respect to the original benchmark problem. 

We consider four different configurations with different channel lengths, see Figure |6j The 
figure also shows the streamlines of the flow, indicating the recirculation. Notice that the non¬ 
symmetry of the streamlines is due to the non-symmetry of the domain. The original length 
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Figure 5. Evolution of the kinetic energy for the ’do-nothing’ and the ’energy’ conditions. 

is chosen according to the ’rule of thumb’. The second computational domain is chosen in 
such a way that the recirculation zone is approximately cut in the middle. Clearly, there is an 
important part of the outflow boundary where the flow is entering the domain. The length of 
the two following domains is increased by 1 meter each. On the longer domain, we have v~ = 0, 
such that our formulation coincides with the ’do-nothing’ condition, but the velocity profile is 
still different from Poiseuille. 





(b) Length 4 m 


(c) Length 5 m 



(d) Length 6 m 


Figure 6 . Tested geometries for the cylinder benchmark with computed streamlines 

In Figure[7]we compare the streamlines of the ’energy’ boundary condition on the small domain 
with the solution on the whole domain. This zoom shows that, although the overall pictures 
look similar, the cut through the vertices clearly has an impact on the flow. Figure [8] shows 
a comparison between the ’energy’ and the ’do nothing’ conditions on the smallest domain. A 
closer look at the re-entrant part of the artificial boundary shows that with the ’energy’ condition, 
the streamlines (lines with circles) are perpendicular to the boundary as expected. This is not 
the case for the streamlines (lines with squares) resulting from the ’do nothing’ condition. 
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Figure 7. Comparison of streamlines for ’energy’ condition: smallest ( circles ) 
and original domain (squares). 



(a) Streamlines 



(b) Zoom 


Figure 8 . Comparison between ’energy’ (circles) and ’do-nothing’ (squares) con¬ 
ditions on the smallest domain 


In order to evaluate the error induced by the boundary condition, we compute the drag coef¬ 
ficient of the cylinder. The reference value, computed on a fine grid with the original geometry, 
is 4.0356. In Table [5j we give the computed values on a sequence of globally refined meshes for 
each of the considered geometries. From the first column, we can see the expected second-order 
convergence starting from hj 4. Already on the coarsest mesh with 160/96/80/64 elements, the 
relative error is below 12% for all geometries. 
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Next we compare the errors due to the boundary conditions on the finest meshes. The relative 
error on the shortest domain is 5.8% and drops to 1.3% and 0.2% on the longer channels. 

Finally, we compare the accuracy of our boundary condition with the ’do-nothing’ condition. 
To this end, Table [5] contains the computed values of the drag coefficients for the latter one. It 
turns out that the relative error for the ’do-nothing’ condition is slightly larger than ours. 



l = 22 

1 = 6 

l = 5 

1 = 4 

l — 4 ’do-nothing’ 

h 

4.492213 

4.493186 

4.476176 

4.374147 

4.370849 

h/2 

3.996568 

4.002665 

3.949923 

3.786517 

3.777821 

h/4 

4.021505 

4.029468 

3.968233 

3.790646 

3.778295 

h/8 

4.032634 

4.041017 

3.978145 

3.797552 

3.784410 

h /16 

4.034720 

4.043192 

3.979928 

3.798635 

3.785252 

hf 32 

4.035179 

4.043660 

3.980291 

3.798886 

3.785401 

Error 

- 

0.2 % 

1.3 % 

5.8 % 

6.2 % 


Table 5. Drag coefficients for different channel length l and mesh resolution h 
and relative error on finest mesh 


5.3. Inviscid flow. In this subsection we consider two test problems for inviscid flows with 
analytical solutions. The first one investigates the numerical dissipation leading to a decrease in 
kinetic energy. The second one is a well-known analytical solution of the Euler equations with 
vortices. 

5.3.1. Standing vortex. The computational domain is D =] — 1, +1[ 2 and the stationary solution 
is given in polar coordinates by v r = 0 and vq — 0 for r > 0.8, vq = 2 — 2.5r for 0.4 < r < 0.8 
and vq — 2.5 r for r < 0.4. Starting with the nodal interpolation of this velocity field as initial 
condition (see Figure [9](a)) , we solve the discrete equations up to t = 5 s on a sequence of uniform 
meshes with number of nodes N = 256,1024,4096,16384 (denoted by meshl,..., mesh4 in the 
following). The numerical dissipation due to the boundary and interior stabilization terms leads 



(a) t = 0s 


(b) t — 5s 


Figure 9. Standing vortex problem on a uniform mesh with N = 1024 nodes: velocity 

to the creation of artificial vortices with very low energy. Since they cannot be detected by the 
vector fields, we use streamlines to visualize this effect in Figure [TO) The temporal behavior of 
the kinetic energy, as well as the numerical dissipation due to the boundary stabilization and 
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Figure 10. Standing vortex problem on a uniform mesh with N = 1024 nodes: streamlines 


interior stabilization terms are shown in Figure [TlJ The loss in kinetic energy on the four meshes 
is: 42%, 11%, 2%, 0.4%. This shows an empirical convergence order between two and three. 

From Figure [Tl|(c) it becomes clear that the SUPG-stabilization is dominant. The ratio be¬ 
tween boundary and interior stabilization at t = 5s for the sequence of meshes is: 25, 221, 876, 
880. 




(a) Kinetic energy. (b) Numerical dissipation (boundary). 



(c) Numerical dissipation (stabilization). 


Figure 11. Standing vortex problem on a sequence of meshes. 
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(b) t « 0.01s 


Figure 12. Comparison between Fraenkeks analytical solution and computation. 

5.3.2. Fraenkel flow. As another example for inviscid flow computation, we consider the analyt¬ 
ical stationary solution of a rotational flow given by Fraenkel in m • This solution with a linear 
far field velocity profile satisfies the wall condition on the lower part of the boundary and presents 
two symmetric vortices before and behind the half-cylinder, see Figure[l2j It has for example been 
used as a test-case in m- We use the computational domain O =] — 3, 3[x]0, 3[\{x^ + x\ < 1}. 



(c) t « 2s 


(d) t w 4s 


Figure 13. Instability of Fraenkel’s analytical solution for inviscid flow: streamlines 


The inflow condition with v D = (^2,0) is given at x\ — —3, the outflow condition is used at 
x\ — 3 and the remaining parts of the boundary are treated as solid walls. Figure [h| shows some 
typical numerical results. Starting from rest, we retrieve the analytical solution after a short 
time t ~ 0.01s. However, starting at £ ^ Is a detachment appears behind the cylinder and the 
solution loses its symmetry and develops other vortices, as one can see in Figure [l3| (c) and (d). 
We have conducted several numerical experiments using different mesh resolutions, larger com¬ 
putational domains, and different time schemes. The results were always quite similar to those 
reported here. We conclude that the analytical solution is not stable. As a further illustration, 
we consider the computation of the flow around a cylinder by reflection of the domain fl. Noting 
that the analytical solution can be prolongated by reflection, we impose the inflow condition 


,,D 


= (|^21,0). As shown in Figure 


14 


the flow loses the symmetry in x\ at about t 
before. The symmetry in X 2 is broken at t ^ 100s. 


Is as 


5.3.3. Jet impact problem. Here we consider the impact of an inviscid jet in order to illustrate 
the influence of the proposed boundary stabilization and the importance of the scaling of the 
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(a) t ~ 0.01s 


(b) t ~ Is 



(c) t ~ 2s (d) t ~ 4s 



Figure 14. Fraenkel flow: solution in the whole domain. 


absolute value matrix. The geometry is a rotated ’T’ of the precise form ft =]0; 2[x]0; 2[x [0; 1] x 
[1; 2] U [0; 1] x [1; 2]. Starting from an inert state, the jet enters with a flat profile (3 — 2 t)t 2 from 
the left, hits a vertical wall and leaves the domain through an upper and a lower channel, as 
shown in Figure |l5(a)[ which is taken at t — 3s; due to symmetry, we only show the solutions 
on half of the domain. 

Comparison with the alternative discretization. We impose the inflow condition on T[ n = Ox ]i;2[, 
the outflow condition on r out =] 1; 2[x 1U] 1; 2[x2 and the wall condition on the remaining parts 
of the boundary r wa n = dQ \ (Ti n U r out ). The flow develops two symmetric large vertices (the 
upper one is shown in Figure |15(a)|), w hich break up into a series of smaller ones and stays 
non-stationary. Figures fl5(b)| and 1 15(c) | show the velocity magnitude computed with strong and 
weak enforcement of the normal velocity, respectively. For the strong implementation we can 
detect oscillations on the upper wall of the boundary, which propagate into the interior of the 
domain. For the weak implementation, oscillations behind the re-entrant corner are also visible, 
but clearly more localized. Notice that both methods employ the same SUPG stabilization, and 
that no shock-capturing terms are used. 
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(a) Flow pattern. (b) Magnitude of v\ (strong). (c) Magnitude of v\ (weak). 

Figure 15. Jet impact problem: strong vs. weak boundary conditions (half domain). 


Influence of the scaling. We use the same test case to investigate the behavior of the discrete 
solutions under scaling. Here, we replace the inflow and outflow boundary conditions by the 
characteristic one with the same v D as before on Tin, v D — 0 on r out and a pressure difference 
equal to 1, i.e. p D = 1 on T[ n and p° = 0 on r out . To this end, we compute the problem in 
the scaled domain Q — sQ with s — 10 and prescribe v D — sv D and p D — s 2 p D , such that the 
continuous solution verifies v — sv and p = s 2 p. The scaling of the domain implies d^ = sd\ l . 
For the sake of clarity, we show in Figure |5~.3. 3 the two computational domains and the pressure 



Figure 16. Jet impact problem: pressure isolines on and Q. 


isolines obtained with the proposed method. In Figures [17(a) and |17(bj| we represent the discrete 
pressures p^ and ph along the symmetry axis y — 1 and y = 10, respectively, at time t — 3. Two 
different scales are used for the sp ace coo rdinates (1 and 10), as well as for the pressure (1 and 
10 2 ), as indicated. The left Figure 17(a) shows the results obtained with the proposed method, 
where 6 is defined in (32). As expected, the invariance to scaling is respected up to machine 
accuracy. For comparison, in the right Figure |17(b)| we show the results obtained with the 
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Figure 17. Jet impact problem: pressure along the symmetry axis. 


alternative discretization (38), which corresponds to 0 = 1 in the definition of the characteristic 
boundary conditions. Clearly, the pressure is no longer invariant under scaling. 


6. Appendix 


We give here the proofs of Lemmas [I] and [2] concerning the spectral decomposition of the Euler 
equations. 


6 . 1 . 


Proof of Lemma 


b 


Let 


x 

y 


be an eigenvector to A. Then 


0x n = Ay, Oyn + (pv n - X)x = 0. 


Multiplying the second equation by On and inserting the first equation gives 

y[0 2 + (pv n - A)A) = 0. 

The case y — 0 immediately gives x — n 1 - in 2D, respectively x = £i, x = £2 in 3D and A = pv n . 
Otherwise we find the two roots A p / m of the quadratic equation 

A 2 — pv n A — 0 2 — 0. 


6.2. Proof of Lemma O We note that 

A p A m = — # 2 , A 2 + A^ = (A p + A m ) 2 — 2A p A m = 2 6 2 + p 2 p 2 
and that \QA n (u)Q\ = R\A\R T . For d — 2, we have 


p\v n \ 0 

0 


4>n ■ n± 

1 f\ A 1 -xA 

0 A p 

0 0 

0 

— A m _ 

, = 

^/02 +A 2 OpK + X) 

1 ( \ rh 4- vA 




so we get 

\^n(u)\e'ip • '0 / = \A\R T Q~ lr i/; • i? T © -1 7// 

= P\ v n\4>n * 0 rt + q 2 + \2'(X + ^P^XA^ + 

_ 5>2 (X + ^m0n)(X / + </>' n )- 
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It goes the same way for d — 3 with (j)^ • q= (j) • <fi f — (frnfin- Using that 


(40) 

we deduce: 


n y"n 

An 


02 + A 2 Q2 + Ap _ Am ^ 4^2 + p 2^2 


URx + A p ^) n )(x / + Ap^) — A, (x + A m ^> n )(x / + A m 0^ 


0 2 + A 2 


1 


A p — A n 


A p — A n 


A p — A n 


1 


\/4fl 2 + p*v. 


P YnJ 0 2 + A2 

(x + Ap<^ n )(x / + A p </4) + (x + A m <^ n )(x / + A m cp' n 
^XX' + (Ap + A m )(0nX / + 4 >' n x) + (Ap + A ^)(p n (p' r 

2XX ; + P v n(<PnX' + 0nX) + + p 2 v n) 4 >n(p'r, 

(X + P v n(pn)(x' + P v n 4 > 'n) + XX 7 + 2 d 2 (p n (p' r 


2„2 
n 


We next compute A n (u) e ip ■ ip' = A R T @ l ip ■ R T @ 1 ip r where 

A“ = 

We obtain: 


P'l-'n 0 0 

0 0 0 

0 0 A m 


A n (u) e 1 p ■ Ip' = pv n (p^ ■ (pf + ^^2 ( A m<N + x)(Am 4>n + X')- 
Finally, noting that 

F 46> 2 + P 2 U 2 ~ 9 + p\v n I, X 2 + (x + PVn(pn ) 2 ^ X 2 + 
we obtain the equivalence of |^4 r i( / ^)|o^ * ^ with 0- 
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